Things that are left to do:

• Handle electrolyte reactive systems (i.e. aqueous systems). The reference system is different so a different activity needs to be used.
• Handle equilibrium systems. I'd imagine a 'cheat' would be to introduce equilibrate in-between steps?
• Create a database of formation energies.

Codecov Report

Attention: Patch coverage is 0% with 155 lines in your changes missing coverage. Please review.

Project coverage is 77.36%. Comparing base (cf96c68) to head (e7b4161).
Report is 7 commits behind head on master.

@@            Coverage Diff             @@
##           master     #284      +/-   ##
==========================================
- Coverage   77.93%   77.36%   -0.58%     
==========================================
  Files         316      321       +5     
  Lines       23585    23796     +211     
==========================================
+ Hits        18382    18410      +28     
- Misses       5203     5386     +183     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Made some headway with the reactive aqueous systems. And got it to work to pure water to get pH:

species = ["water","hydronium","hydroxide"]
reactions = [("water"=>-2,"hydronium"=>1,"hydroxide"=>1)]

eosmodel = SAFTVREMie(["water"],["hydronium","hydroxide"])
model = ReactiveAqModel(species,reactions; model = eosmodel)

Clapeyron.pH(model,1e5, 298.15, [1.0, 0.0, 0.0]; z0=[1e-9]) # 7.0978252141976785

The difficult part here is the reference systems...

In the case of electrolytes, since the ions can't exist in isolation, the reference formation energies are measured at 1 mol/kg of the salt. This is fine when there is a single reaction (2 H2O<-> H3O+ + OH-) since that would just be the H2O 'salt'. However, if you have a weak acid, and thus two reactions (HA + H2O <-> H3O+ + A-), there are three ions, H3O+, A- and OH-. What's the reference system in this case? My feeling is I'll treat the reference as reaction-dependent. So the reference will change for each reaction, but at least the way the reference is treated will be consistent.

Maybe we could differentiate references by parameter name? (Gf_aq fvs Hf0) there are a finite amount of references anyway. I suppose we can convert between those if we have info of at least one of them

So at the moment I'm differentiating them by the model type (ReactiveModel vs ReactiveAqModel). I think since the way the equations are solved and the concentration scales involved are different, it's best to keep it that way. However, I seem to have failed at defining the struct types since, ideally, ReactiveAqModel <: ReactiveModel, but I couldn't get it to work. Any advice @longemen3000?

i was just looking at the code, and the main difference between the aqueous and normal reactive model are the value of μ_ref. For example, activity_coefficient and activity can be defined as a function of μ_mix - μ_ref . i added an API that allows this explicit separation, maybe we could define transformations of references by modifying the value of μ_ref (at least for pure -> aqueous)